home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / mmid.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  1.3 KB  |  57 lines
  1. /*
  2. ### modified mid-point method ###
  3. */
  4.  
  5. extern int kount;
  6. void mmid(y,dydx,nvar,xs,htot,nstep,yout,int_option)
  7. double y[],dydx[],xs,htot,yout[];
  8. int nvar, nstep,int_option;
  9. {
  10.     int n,i;
  11.     double x,swap,h2,h,*yt,*ym,*yn,*vfull,*dvector();
  12.     extern int model,polar_coord,full_dim;
  13.     extern double *param;
  14.     extern int (*f_p)();
  15.  
  16.  
  17.     yt=dvector(0,nvar-1);
  18.     ym=dvector(0,nvar-1);
  19.     yn=dvector(0,nvar-1);
  20.     vfull=dvector(0,full_dim-1);
  21.     h = htot/nstep;
  22.     for(i=0;i<nvar;i++){
  23.         ym[i] = y[i];
  24.         yn[i] = y[i] + h * dydx[i];
  25.     }
  26.     x = xs + h;
  27.     (int) f_p(yout,0,yn,param,x,nvar);
  28.     if(int_option==1) {
  29.         for (i=0;i<nvar;i++)
  30.             yt[i] = 0.5 * (ym[i] + yn[i] + h * yout[i]);
  31.         to_full_variables(vfull,yt,polar_coord);
  32.         (void) draw_record_orbit(vfull,kount++,0);
  33.     }
  34.     h2 = 2.0 * h;
  35.     for(n=1;n<nstep;n++){
  36.         for(i=0;i<nvar;i++){
  37.             swap = ym[i] + h2 * yout[i];
  38.             ym[i] = yn[i];
  39.             yn[i] = swap;
  40.         }
  41.         x += h;
  42.         (int) f_p(yout,0,yn,param,x,nvar);
  43.         if(int_option==1) {
  44.             for (i=0;i<nvar;i++)
  45.                 yt[i] = 0.5 * (ym[i] + yn[i] + h * yout[i]);
  46.             to_full_variables(vfull,yt,polar_coord);
  47.             (void) draw_record_orbit(vfull,kount++,0);
  48.         }
  49.     }
  50.     for (i=0;i<nvar;i++)
  51.         yout[i] = 0.5 * (ym[i] + yn[i] + h * yout[i]);
  52.     (void) free_dvector(vfull,0,full_dim-1);
  53.     (void) free_dvector(yn,0,nvar-1);
  54.     (void) free_dvector(ym,0,nvar-1);
  55.     (void) free_dvector(yt,0,nvar-1);
  56. }
  57.